function ddelta = rcnl_ddelta( ...
        Pi,           ...
        rho,          ...
        delta,        ...
        xnonlin,      ...
        cdindex,      ...
        dfull,        ...
        rho_transform ...
    )

    mu         = (xnonlin * Pi) .* dfull;
    cdindex_   = [0; cdindex];
    oneInvRho  = 1 / (1 - rho);
    if rho_transform
        dtransform = rho * (1 - rho/rho_transform);
    else
        dtransform = 1;
    end

    ddelta = zeros([size(xnonlin, 1) size(xnonlin, 2) + 1]);
    for m = 1:length(cdindex)
        l = cdindex_(m)+1;
        r = cdindex_(m+1);

        delta_m   = delta(l:r);
        mu_m      = mu(l:r, :);
        xnonlin_m = xnonlin(l:r, :);
        dfull_m   = dfull(l:r, :);

        num       = exp(oneInvRho * (delta_m + mu_m));
        den       = sum(num, 1);
        s_ijrt_g  = num ./ den;
        I_irt_g   = den.^(1 - rho);
        s_irt_g   = (I_irt_g ./ (1 + I_irt_g));
        s_ijrt    = s_ijrt_g .* s_irt_g;
        s_jrt     = mean(s_ijrt, 2);
        N_rt      = size(mu_m, 2);

        D_delta = oneInvRho * (diag(s_jrt) - s_ijrt * (rho * s_ijrt_g + (1 - rho) * s_ijrt)' / N_rt);

        i_aux = 1:size(xnonlin_m, 2);
        D_aux = (oneInvRho * dfull_m .* s_ijrt);
        X_aux = xnonlin_m' * (rho * s_ijrt_g + (1 - rho) * s_ijrt);
        D_nonlin = splitapply(@(i) mean(D_aux .* (xnonlin_m(:, i) - X_aux(i, :)), 2), i_aux, i_aux);

        D_rho = oneInvRho * oneInvRho * mean( ...
            s_ijrt .* ( ...
                (delta_m + mu_m) ...
                - sum((delta_m + mu_m) .* (rho * s_ijrt_g + (1 - rho) * s_ijrt), 1) ...
                - (1 - rho) * (1 - rho) * log(den) .* (1 - s_irt_g) ...
            ), 2 ...
        );

        ddelta(l:r, :) = - D_delta \ [D_nonlin D_rho * dtransform];
    end
end
